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SUMMARY 


The interaction of a planar shock wave with one or more vortexes 
is computed using a pseudospectral method and a finite difference 
method. The paper emphasizes the development of the spectral method. 

In both methods the shock wave is fitted as a boundary of the computa- 
tional domain. The results show very good agreement between ^oth 
computational methods. The spectral method is, however, restricted to 
smaller time steps and requires use of filtering techniques. 

INTRODUCTION 

In their recent paper, Pao and Salas (ref. 1) presented a finite 
difference solution to the Euler equations governing the phenomena of 
shock wave interaction with an isolated vortex. Their study emphasized 
the acoustic aspects of the problem. Zang, Hussaini, and Bushnell 
(ref. 2) extended this numerical approach to the problem of shock wave 
interaction with entropy spottiness (or hot spots), plane waves, etc. 
The present study is a continuation of this effort which concentrates 
on gaining Insight into the nonlinear dynamics of the transient pro- 
cesses involved in the passage of a shock wave over a single vortex, a 
vortex street, and a hot spot. The governing equations are solved by 
two different numerical methods j the well known second— order j finite 
difference method originated by MacCormack, and a pseudospectral method 




of high accuracy. The finite difference solutions are calculated on a 
very fine grid and are used for comparison with the solutions obtained 
with the spectral method. The main emphasis of the paper is on 
developing a viable spectral technique. 

Spectral methods have been demonstrated (refs. 3,4) to be powerful 
alternatives to finite difference methods for the numerical solution of 
smooth flows. Recently, the works of Gottlieb, Lustman, and Orszag 
(ref. 5) and Zang and Hussaini (ref. 6) have shown their applicability 
to one-dimensional compressible flows with shocks. However, for multi- 
dimensional flows in the presence of shock waves, spectral methods have 
not been developed to the same level of proficiency as finite difference 
methods. The present paper is thus an attempt to break new ground in 
this area. 

STATEMENT OF THE PROBLEM 

The physical problem that we model, the shock-vortex interaction 
problem, corresponds to an infinite, initially planar normal shock wave 
moving from left to right into a downstream region containing a flow 
field representative of one or more vortexes, or a hot spot. Since the 
shock wave travels at a speed greater than the local downstream speed 
of sound, the flow downstream of the shock wave remains unaware of the 
oncoming shock until it is overtaken by it. In order to model the 
interaction of the shock wave with some given flow field ahead of it, 
it is only necessary to compute the flow field upstream of the shock. 

The computational domain, therefore, consists of the region between 
some left boundary, judiciously chosen such that it will be far from 
the interaction region, and the shock wave front itself. A coordinate 
transformation is used to map the infinite lateral extent onto a finite 
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domain, thus concentrating the available grid points in the vicinity of 
the interaction. 

If the relative shock Mach number, M^, is sufficiently high 

(M > 2.08), the flow upstream of the shock remains supersonic and 

s 

signals from the shock-vortex interaction will not reach the left boun- 
dary. In this case, the left boundary corresponds to a supersonic 
inflow, and all dependent variables can be prescribed on it. Most of 
the cases to be studied here will be of this type. However, if the 
relative shock Mach number is low, then the first signal from the 
interaction will reach the left boundary at some finite time t^. For 
times greater than t^, special procedures are required to prevent 
contamination of the solution in the region of interest by spurious 
wave reflections from the left boundary. Two approaches are available 
to achieve this end. Either radiation type boundary conditions can be 
applied at the left boundary, or the boundary can be allowed to move in 
such a way as to insure that it is inaccessible to the initial signal 
emanating from the interaction region. Both approaches have been tried 
for this problem and in a comparison with the experiment of reference 7 
we will use the first method for illustration. 

On the right, the computational region is bounded by the shock 
wave. Downstream of the shock the flow field is given analytically. 

The flow field immediately upstream of the shock, as well as the shape 
and velocity of the shock, are evaluated such that the Rankine-Hugoniot 
Jump conditions and the compatibility condition reaching the shock wave 
from the upstream side are simultaneously satisfied. 
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GOVERNING EQUATIONS 


The unsteady, two-dimensional, compressible, Euler equations are 
written in the form. 


Qt + A’Q^ + B’Qy = 0 


where 


Q = [^P,u,v,S 2 
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with subscripts denoting partial derivatives with respect to the inde- 
pendent variables. The natural logarithm of the pressure, the speed of 
sound, and the entropy are represented by P, a, and S, respectively, 
and Y is the isentropic exponent. The velocity in the x and y direc- 
tions are u and v, respectively. All variables are normalized with 
respect to reference conditions at downstream infinity, as in reference 

1 . 

The physical domain is defined by 
h(t) < X < Xg(y,t) 

( 6 ) 


— 00 ^ y ^ 00 
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where h and x are the abscissas of the left boundary and the shook. 

s 

respectively . 

The physical domain is mapped onto a computational domain by the 
following transformation, 


X - h(t) _ 

(7) 

^ Xg(y,t) - h(t) 


tanh(ay) + 1 
^ " 2 

(8) 

T = t 

(9) 


The computational domain is thus 

0 ^ X ^ 1 
0 ^ Y ^ 1 

and the governing equations (1-4) are rewritten as 
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here U and V are the contravariant velocity components defined by 


U = uX^ + vXy 


(15) 


V = + uY^ + vYy 


(16) 


SOLUTION TECHNIQUE 

Let k denote the time level and let At be the time step incre- 
ment, the time discretization of equation (12) is then as follows. 


Q* 


= [l - At (L^ + Ly)]q‘ 

= |[q^ + (1 - At (L^ + Ly))Q*] 


(17) 

( 18 ) 


where the spatial operators L^^, and L^ represent approximations to 

A8 and B8y, respectively. In the finite difference MacCormack method, 
“X ^ 

these operators are evaluated as two points forward and two points back- 
ward differences in the predictor (Eq. (17)) and corrector (Eq. (18)) 
levels, respectively. In the pseudospectral method studied here, 
the solution Q is first expanded as a double Chebyshev series, 

p=0 q=0 
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Q(X,Y,T) = 



where 


5 = 2X - 1 (20) 

n = 2Y - 1 (21) 

and T and t are the Chebyshev polynomials of degree p and q 
P q 

defined by 

Tp(0 = cos(p cos^(5)) (22) 

X (n) = cos(q cos^(Ti)) (23) 

Q. 

The (M+1).(N+1) Chebyshev coefficients Qp^ are such that at 

(Vn'T)' 

•M N 
p=0 q=0 

and 

= I = I 1) '2®> 

' I ^ 1) = I '2®> 

The derivatives appearing in the operators and Ly are then 
evaluated as 

«X= 2 L X %q ^ t, (27) 

p=0 q=0 ^ ’ 

The last relation can be expressed by another Chebyshev series 
M N , 

«x= 2 X;X 9pq "p 7, (28) 

p=0 q=0 
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where 


M 


and 


Q = 

Tq 


P m-ptl 


“V 


m+p odd 


Cq 2 

c = 1, P > 0 
P 


(29) 


(30) 

(31) 


The Qy derivative is evaluated in a similar fashion. 

Assuming the characteristic wave speed to be X, the time step 

for the finite difference method with uniformally spaced N grid 

points over a unit interval is At < while for the spectral 

method with N modes the time step is At < This strong restric- 

N X 

tion on the time step is at present one of the main disadvantages of 
the spectral approach. Another disadvantage is its tendency to 
develop spurious oscillations. The exact cause of these oscillations 
is presentdy unknown, but its instantaneous spreading over the field 
must be caused by the global nature of the spectral method. Various 
filtering techniques are, however, available in order to recover the 
solution. The results presented here were obtained by applying a 
von Hann Window filter, see reference 6 for details, at every 160 
time steps. A comprehensive study of various filters, as well as 
of the effect of filtering on the solution is under preparation. 

The evaluation of the shock wave shape and velocity followed 
the same procedure described in Ref. 1, except that in the spectral 
formulation, the derivatives that must be evaluated on the downstream 
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side of the shock are expressed as Chebyshev expansions as given 
by Eq. (28). 

At the left boundary, all variables were specified for supersonic 
inflow. For the case of subsonic inflow, the two velocity components 
and the entropy were specified, while the pressure was computed from 
a quasi-one-dimensional characteristic as described in Ref. 8. 

RESULTS 

Perhaps the simplest interaction to consider is that of a planar 
shock wave with a hot spot, as shown in figure 1. The flow field 
downstream of the shock wave situated at x=0 at t=0 is taken as a 
quiescent field whose temperature, a, distribution is given by 

a = exp (-((x-Xq)^ + (y-y^)^ )/ 2r^) (32) 

where for this case k = 0.25, r = .125, x^ = 0.5, y^ = 0.0, and 

M = 3 at t = 0. The flow field and vorticity distribution obtained 
s 

by the finite difference method and the spectral method after the 
shock wave has passed over the hot spot is shown in figure 2. The 
finite difference solution presented here, and in all cases that 
follow, was obtained with 75 mesh; points in the X direction and 60 
mesh points in the Y direction. The spectral solutions were all 
obtained with 32 collocation points in the X direction and 16 colloca- 
tion points in the Y direction. Qualitatively and quantitatively, 
there is very little difference between the two solutions. Both show 
the two counterrotating vortexes downstream of the shock, typical of 
this interaction. 

Figure 3 shows the velocity field for a single vortex about to 
interact with a shock wave of the same initial strength as in the 
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previous case. The upstream conditions here are obtained by assuming 
a constant density field, calculating the velocity from the stream 
function, 

logyr^+(x-x^)^+ (y-y^)^, (33) 

the pressure from Bernoulli's relation, and the temperature from , 
the equation of state. For the case shown in figure 3, the circula- 
tion k= 2 and the softening scale r=0.1. This model approaches 
an idealized incompressible point vortex at large distances from its 
center but is much smoother near the center. Figure 4 shows the 
resulting pressure field at two different time steps. At the later 
time the shock wave has passed over the vortex. For this case, the 
spectral method seems to be resolving the pressure field more accura- 
tely than the finite difference method. Overall, the results are 
qualitatively very similar. Finally, figures 5 and 6 show the 
results for the interaction with a Karman vortex street simulating 
the conditions of the experiment reported in Ref. 7. For this case, 
the stream function representing the vortex field is given by 

^ log(A/B) (34) 

where 

A = cosh + (y - ^b)^ ) - cos (^ (x - |-c)) (35) 

B = cosh (^^r^ t (y + ^b)^) - cos (^ (x + |c)) (36) 

To match the experiment, the circulation, core radius, shock Mach 
number and vortex separation parameters were determined as k= 0.186, 
r= 0.1, Mg = 1.3, c= .33, and b = 0.17. For this calculation, the 


10 



inflow Mach number was subsonic and radiation boundary conditions were 
applied at the left boundary. The results shown in figure 6 are in 
agreement with the experimentally observed (ref. 7) longitudinal com- 
pression and lateral elongation of the vortex field after passage 
through the shock. 

CONCLUSION 

A pseudospectral technique has been presented for the shock-vortex 
interaction problem. The solutions obtained are in good agreement with 
those predicted by a standard finite difference method. The proper 
treatment of the shock wave, which is one of the major problems with 
any numerical technique, was resolved by fitting the shock wave as a 
boundary of the computational plane. The spectral technique seems 
capable of resolving features of the flow field which are difficult to 
resolve by finite difference. The spectral method is, however, 
presently restricted to very small time steps, and it requires the 
use of filtering techniques which need further investigation. 
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Figure 1.- Surface plot of temperature distribution for a hot 

spot about to be overtaken by a shock wave moving 

at = 3. 
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Figure 3.- Velocity field depicting single vortex about to 
interact with a shock wave. 
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